const PI = 3.1415926535898;
const a = 6378245;  //长半径
const f = 298.3;   //扁率的倒数 （扁率：(a-b)/a）;
const e2 = 1 - ((f - 1) / f) * ((f - 1) / f);        //第一偏心率的平方
const e12 = (f / (f - 1)) * (f / (f - 1)) - 1;       //第二偏心率的平方
let L0;
/**
 * 角度转弧度
 */
function degree2Rad(degree) {
    let Sign, Rad;
    if (degree >= 0) {
        Sign = 1;
    }
    else {
        Sign = -1;
    }
    degree = Math.abs(degree); //绝对值
    Rad = Sign * degree * PI / 180.0;
    return Rad;
}

/**
 * 弧度转角度
 * @param rad
 * @return {number}
 */
function rad2Degree(rad) {
    let degree, Sign;
    if (rad >= 0) {
        Sign = 1;
    }
    else {
        Sign = -1;
    }
    rad = Math.abs(rad * 180.0 / PI);
    degree = Sign * rad;
    return degree;
}

/**
 * WGS84 转 北京54 坐标系
 * @param lon
 * @param lat
 * @return {{x: number, y: number}}
 * @constructor
 */
function WTB(lon, lat) {
    let MyL0;  //中央子午线
    let A1, A2, A3, A4;
    let B, L;  //中央子午线弧度、纬度弧度、经度弧度
    let X;   //由赤道至纬度为B的子午线弧长   （P106   5-41）
    let N;   //椭球的卯酉圈曲率半径
    let t, t2, m, m2, ng2, cosB, sinB;

    // 克拉索夫斯基椭球
    A1 = 111134.8611;
    A2 = -16036.4803;
    A3 = 16.8281;
    A4 = -0.0220;

    //计算当地中央子午线，3度带
    MyL0 = Math.abs(lon);
    MyL0 = Math.floor(MyL0);
    MyL0 = 6 * Math.floor(MyL0 / 6) - 3;
    L0 = degree2Rad(MyL0); //将中央子午线转换为弧度
    B = degree2Rad(lat);
    L = degree2Rad(lon);

    X = A1 * B * 180.0 / PI + A2 * Math.sin(2 * B) + A3 * Math.sin(4 * B) + A4 * Math.sin(6 * B);
    sinB = Math.sin(B);
    cosB = Math.cos(B);
    t = Math.tan(B);
    t2 = t * t;

    N = a / Math.sqrt(1 - e2 * sinB * sinB);
    m = cosB * (L - L0);
    m2 = m * m;
    ng2 = cosB * cosB * e2 / (1 - e2);
    let x = X + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2);
    let y = N * m * (1 + m2 * ((1 - t2 + ng2) / 6.0 + m2 * (5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2) / 120.0));

    return {
        x: x,
        y: y
    };
}

/**
 * 北京54 转 WGS84 坐标系
 * @param x
 * @param y
 * @return {{lon: number, lat: number}}
 * @constructor
 */
function BTW(x, y) {

    let N;   //椭球的卯酉圈曲率半径
    let t, t2, ng2, cosB, sinB, V, yN, preB0, B0, eta;
    let A1, A2, A3, A4;
    let sgp_lat, sgp_lon;

    // 克拉索夫斯基椭球
    A1 = 111134.8611;
    A2 = -16036.4803;
    A3 = 16.8281;
    A4 = -0.0220;

    B0 = x / A1;
    do {
        preB0 = B0;
        B0 = B0 * PI / 180.0;
        B0 = (x - (A2 * Math.sin(2 * B0) + A3 * Math.sin(4 * B0) + A4 * Math.sin(6 * B0))) / A1;
        eta = Math.abs(B0 - preB0);
    } while (eta > 0.000000001);

    B0 = B0 * PI / 180.0;
    sinB = Math.sin(B0);
    cosB = Math.cos(B0);
    t = Math.tan(B0);
    t2 = t * t;
    N = a / Math.sqrt(1 - e2 * sinB * sinB);
    ng2 = cosB * cosB * e2 / (1 - e2);
    V = Math.sqrt(1 + ng2);
    yN = y / N;

    sgp_lat = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2;
    sgp_lon = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB;

    return {
        lon: rad2Degree(sgp_lon),
        lat: rad2Degree(sgp_lat)
    }
}

/**
 * 高斯转经纬度
 * @param X
 * @param Y
 * @return {{latitude: Number, longitude: Number}}
 * @constructor
 */
function GaussToBL(X, Y) {
    let ProjNo;
    let ZoneWide; //带宽
    let output = new Array(2);
    let longitude1, latitude1, longitude0, X0, Y0, xval, yval;//latitude0,
    let e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
    iPI = 3.14159265358979324 / 180.0; //3.1415926535898/180.0;
    // a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数
    a = 6378140.0; f = 1.0 / 298.257; //80年西安坐标系参数
    ZoneWide = 6; //6度带宽
    ProjNo = parseInt(X / 1000000); //查找带号
    longitude0 = (ProjNo - 1) * ZoneWide + ZoneWide / 2;
    longitude0 = longitude0 * iPI; //中央经线


    X0 = ProjNo * 1000000 + 500000;
    Y0 = 0;
    xval = X - X0; yval = Y - Y0; //带内大地坐标
    e2 = 2 * f - f * f;
    e1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2));
    ee = e2 / (1 - e2);
    M = yval;
    u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
    fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.sin(4 * u) + (151 * e1 * e1 * e1 / 96) *
        Math.sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.sin(8 * u);
    C = ee * Math.cos(fai) * Math.cos(fai);
    T = Math.tan(fai) * Math.tan(fai);
    NN = a / Math.sqrt(1.0 - e2 * Math.sin(fai) * Math.sin(fai));
    R = a * (1 - e2) / Math.sqrt((1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin
        (fai) * Math.sin(fai)));
    D = xval / NN;
    //计算经度(Longitude) 纬度(Latitude)
    longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
        * D * D * D * D / 120) / Math.cos(fai);
    latitude1 = fai - (NN * Math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
        + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
    //转换为度 DD
    output[0] = longitude1 / iPI;
    output[1] = latitude1 / iPI;
    return {
        longitude: output[0],
        latitude: output[1]
    };
}

/**
 * 经纬度转高斯
 * @param longitude
 * @param latitude
 * @return {{x: Number, y: Number}}
 * @constructor
 */
function BLToGauss(longitude, latitude) {
    let ProjNo = 0;
    let ZoneWide; //带宽
    let ret = Array(2);
    let longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
    let a, f, e2, ee, NN, T, C, A, M, iPI;
    iPI = 0.0174532925199433;// 3.1415926535898/180.0;
    ZoneWide = 6;// 6度带宽
    // a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
    a = 6378140.0; f = 1 / 298.257; //80年西安坐标系参数
    ProjNo = parseInt(longitude / ZoneWide);
    longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
    longitude0 = longitude0 * iPI;
    latitude0 = 0;
    longitude1 = longitude * iPI; //经度转换为弧度
    latitude1 = latitude * iPI; //纬度转换为弧度
    e2 = 2 * f - f * f;
    ee = e2 * (1.0 - e2);
    NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1) * Math.sin(latitude1));
    T = Math.tan(latitude1) * Math.tan(latitude1);
    C = ee * Math.cos(latitude1) * Math.cos(latitude1);
    A = (longitude1 - longitude0) * Math.cos(latitude1);
    M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2
        * e2 / 1024) * Math.sin(2 * latitude1)
        + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude1));
    xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
    yval = M + NN * Math.tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
        + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
    X0 = 1000000 * (ProjNo + 1) + 500000;
    Y0 = 0;
    xval = xval + X0; yval = yval + Y0;
    ret[0] = xval;
    ret[1] = yval;

    return {
        x: ret[0],
        y: ret[1]
    };
}


export default {
    WGS84ToBeiJin54(lon, lat) {
        return WTB(lon, lat);
    },
    BeiJin54ToWGS84(x, y) {
        return BTW(x, y);
    },
    LatLonToGauss(lon, lat){
        return BLToGauss(lon, lat);
    },
    GaussToLonLat(X, Y) {
        return GaussToBL(X, Y);
    }
}
